// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//  File: avtSPHResampleFilter.C
// ************************************************************************* //
#include <algorithm>

#include <vtkDataSet.h>
#include <avtSPHResampleFilter.h>
#include <avtParallel.h>
#include <avtResampleSelection.h>

#include <vtkCellData.h>
#include <vtkCellType.h>
#include <vtkPoints.h>
#include <vtkPointData.h>
#include <vtkDoubleArray.h>
#include <vtkRectilinearGrid.h>
#include <vtkAppendPolyData.h>
#include <vtkCleanPolyData.h>
#include <vtkPolygon.h>
#include <vtkMath.h>

#include <InvalidVariableException.h>
#include <DebugStream.h>
#include <TimingsManager.h>

#ifdef PARALLEL
#include <mpi.h>
#endif

using std::vector;
using std::min;
using std::max;

const double KERNEL_EXTENT = 2.0;

// ****************************************************************************
//  Method: avtSPHResampleFilter constructor
//
//  Programmer: harrison37 -- generated by xml2avt
//  Creation:   Fri Dec 5 13:50:31 PST 2014
//
// ****************************************************************************

avtSPHResampleFilter::avtSPHResampleFilter()
{
    keepNodeZone = false;
    nDim = 0;
}


// ****************************************************************************
//  Method: avtSPHResampleFilter destructor
//
//  Programmer: harrison37 -- generated by xml2avt
//  Creation:   Fri Dec 5 13:50:31 PST 2014
//
//  Modifications:
//
// ****************************************************************************

avtSPHResampleFilter::~avtSPHResampleFilter()
{
}


// ****************************************************************************
//  Method:  avtSPHResampleFilter::Create
//
//  Programmer: harrison37 -- generated by xml2avt
//  Creation:   Fri Dec 5 13:50:31 PST 2014
//
// ****************************************************************************

avtFilter *
avtSPHResampleFilter::Create()
{
    return new avtSPHResampleFilter();
}


// ****************************************************************************
//  Method:      avtSPHResampleFilter::SetAtts
//
//  Purpose:
//      Sets the state of the filter based on the attribute object.
//
//  Arguments:
//      a        The attributes to use.
//
//  Programmer: harrison37 -- generated by xml2avt
//  Creation:   Fri Dec 5 13:50:31 PST 2014
//
// ****************************************************************************

void
avtSPHResampleFilter::SetAtts(const AttributeGroup *a)
{
    atts = *(const SPHResampleAttributes*)a;
}


// ****************************************************************************
//  Method: avtSPHResampleFilter::Equivalent
//
//  Purpose:
//      Returns true if creating a new avtSPHResampleFilter with the given
//      parameters would result in an equivalent avtSPHResampleFilter.
//
//  Programmer: harrison37 -- generated by xml2avt
//  Creation:   Fri Dec 5 13:50:31 PST 2014
//
// ****************************************************************************

bool
avtSPHResampleFilter::Equivalent(const AttributeGroup *a)
{
    return (atts == *(SPHResampleAttributes*)a);
}

// ****************************************************************************
// Method:  avtSPHResampleFilter::kernelValue
//
// Purpose:
//   Return the value of the kernel for a given normalized distance.
//
// Arguments:
//   etaMagnitude   The normalized distance (r/h)
//   Hdet           Determinant of the support tensor (if the tensor is symmetric,
//                  this is equivalent to 1/h^3 and can never be zeror under
//                  normal circumstances)
//
// Programmer:  Cody Raskin (raskin1@llnl.gov)
// Creation:    July 1, 2015
//
// ****************************************************************************

template<int Dim>
double
avtSPHResampleFilter::kernelValue(double etaMagnitude, double Hdet) {
    double result = 0;
    double volumeNormalization = 10.0/(7.0*3.14159);
    if(Dim==3)
        volumeNormalization = 1.0/3.14159;
    
    if (etaMagnitude < 1.0) {
        double eta2 = etaMagnitude*etaMagnitude;
        result = volumeNormalization*Hdet*(1.0 - 1.5*eta2 + 0.75*eta2*etaMagnitude);
    } else if (etaMagnitude < 2.0) {
        result = volumeNormalization*Hdet*0.25*pow(2.0 - etaMagnitude,3.0);
    } else {
        result = 0;
    }
    
    return result;
}

// ****************************************************************************
// Method:  avtSPHResampleFilter::kernelGradValue
//
// Purpose:
//   Return the value of the gradient of the kernel for a given normalized distance
//
// Arguments:
//   etaMagnitude   The normalized distance (r/h)
//   Hdet           Determinant of the support tensor (if the tensor is symmetric,
//                  this is equivalent to 1/h^3 and can never be zero under normal
//                  circumstances)
//
// Programmer:  Cody Raskin
// Creation:    July 1, 2015
//
// ****************************************************************************

template<int Dim>
double
avtSPHResampleFilter::kernelGradValue(double etaMagnitude, double Hdet) {
    
    double volumeNormalization = 10.0/(7.0*3.14159);
    if(Dim==3)
        volumeNormalization = 1.0/3.14159;
    
    if (etaMagnitude < 1.0) {
        return -volumeNormalization*Hdet*(3.0 - 2.25*etaMagnitude)*etaMagnitude;
    } else if (etaMagnitude < 2.0) {
        return -volumeNormalization*Hdet*0.75*pow(2.0 - etaMagnitude,2.0);
    } else {
        return 0.0;
    }
}

// ****************************************************************************
// Method:  stepSize
//
// Purpose:
//   Compute the step size of the grid.
//
// Arguments:
//   minXYZ       Vector for the min grid location in each direction
//   maxXYZ       Vector for the max grid location in each direction
//   nsample    Vector for the number of grid cells in each direction
//
// Programmer:  Cody Raskin
// Creation:    July 1, 2015
//
// ****************************************************************************

template<int Dim>
inline
vector<double>
stepSize(vector<double>& minXYZ, vector<double>& maxXYZ, vector<int>& nsample) {
    vector<double> result(Dim);
    
    for (int i = 0; i < Dim; ++i) {
        result[i] = (maxXYZ[i]-minXYZ[i]) / nsample[i];
    }
    
    return result;
}

vector<double>
avtSPHResampleFilter::GetStepSizeXYZ()
{
    // Global Lattice Info to calculate step sizes
    vector<double> globalLatticeMin(nDim);
    vector<double> globalLatticeMax(nDim);
    vector<int> globalLatticeSamples(nDim);
    
    globalLatticeSamples[0] = atts.GetXnum();
    globalLatticeSamples[1] = atts.GetYnum();
    globalLatticeMin[0] = atts.GetMinX();
    globalLatticeMin[1] = atts.GetMinY();
    globalLatticeMax[0] = atts.GetMaxX();
    globalLatticeMax[1] = atts.GetMaxY();
    
    if(nDim == 3)
    {
        globalLatticeMin[2] = atts.GetMinZ();
        globalLatticeMax[2] = atts.GetMaxZ();
        globalLatticeSamples[2] = atts.GetZnum();
    }
    
    vector<double> stepXYZ = (nDim == 3) ? stepSize<3>(globalLatticeMin, globalLatticeMax, globalLatticeSamples)
                                         : stepSize<2>(globalLatticeMin, globalLatticeMax, globalLatticeSamples);
    
    return stepXYZ;
}

// ****************************************************************************
// Method:  makeIndex
//
// Purpose:
//   Generate the dimension appropriate index vector from integers.
//
// Arguments:
//   ix Dimension 1 index
//   iy Dimension 2 index
//   iz Dimension 3 index
//
// Programmer:  Cody Raskin
// Creation:    July 1, 2015
//
// ****************************************************************************

template<int Dim>
inline
vector<int>
makeIndex(const int ix,
          const int iy,
          const int iz) {
    //REQUIRE(false);
}

template<>
inline
vector<int>
makeIndex<1>(const int ix,
             const int iy,
             const int iz) {
    vector<int> result(1);
    result[0] = ix;
    return result;
}

template<>
inline
vector<int>
makeIndex<2>(const int ix,
             const int iy,
             const int iz) {
    vector<int> result(2);
    result[0] = ix;
    result[1] = iy;
    return result;
}

template<>
inline
vector<int>
makeIndex<3>(const int ix,
             const int iy,
             const int iz) {
    vector<int> result(3);
    result[0] = ix;
    result[1] = iy;
    result[2] = iz;
    return result;
}

// ****************************************************************************
// Method:  latticePoints
//
// Purpose:
//   Compute the indicies of the lattice over which a given particle overlaps.
//
// Arguments:
//   ri             Location of the particle
//   extent         Extent of the kernel
//   minXYZ         Vector for the min grid location in each direction
//   stepSizeXYZ    Vector describing the distance between grid location in each direction
//   gridCellsXYZ   Vector for the number of grid cells in each direction
//
// Programmer:  Cody Raskin
// Creation:    July 1, 2015
//
// ****************************************************************************

template<int Dim>
inline
vector< vector<int> >
latticePoints(const double* ri,
              vector<double>& extent,
              vector<double>& minXYZ,
              vector<double>& stepSizeXYZ,
              const vector<int>& gridCellsXYZ) {
    
    // Find the min & max indicies in each dimension.
    vector<int> imin(Dim);
    vector<int> imax(Dim);
    int ntot = 1;
    
    for (int i = 0; i < Dim; ++i) {
        imin[i] = max(0, min(gridCellsXYZ[i] - 1, int((ri[i] - extent[i] - minXYZ[i])/stepSizeXYZ[i])));
        imax[i] = max(0, min(gridCellsXYZ[i] - 1, int((ri[i] + extent[i] - minXYZ[i])/stepSizeXYZ[i]) + 1));
        ntot *= imax[i] - imin[i] + 1;
    }
    
    // Now build the result with the full set of indicies.
    vector< vector<int> > result;
    result.reserve(ntot);
    int ix = imin[0];
    int iy = imin[1];
    int iz = Dim == 3 ? imin[2] : 0;
    int i = 0;
    
    while (i < ntot) {
        result.push_back(makeIndex<Dim>(ix, iy, iz));
        ++i;
        ++ix;
        if (ix > imax[0]) {
            ix = imin[0];
            ++iy;
            if (iy > imax[1]) {
                iy = imin[1];
                ++iz;
            }
        }
    }
    
    return result;
}

// ****************************************************************************
// Method:  latticePosition
//
// Purpose:
//   Compute the position for the given lattice index.
//
// Arguments:
//   cellIndexXYZ   Vector of grid indices
//   minXYZ         Vector for the minimum grid location in each direction
//   stepSizeXYZ    Vector describing the distance between grid locations in
//                  each direction
//
// Programmer:  Cody Raskin
// Creation:    July 1, 2015
//
// Modifications:
//
//      Kevin Griffin, Wed May 25 14:38:01 PDT 2016
//      Changed the addition of 0.5*... to the subtraction of 0.5*... to get
//      the correct cell center position. Renamed the method parameters and
//      removed the unused latticeMax parameter.
//
// ****************************************************************************

template<int Dim>
inline
vector<double>
latticePosition(const vector<int>& cellIndexXYZ, vector<double>& minXYZ, vector<double>& stepSizeXYZ)
{
    vector<double> result(Dim);
    for(int k=0;k<Dim;++k) {
        result[k] = minXYZ[k]+cellIndexXYZ[k]*stepSizeXYZ[k] - 0.5*stepSizeXYZ[k];
    }
    
    return result;
}

// ****************************************************************************
// Method:  latticeIndex
//
// Purpose:
//    Convert the tuple of indicies into the flat index for the lattice arrays.
//
// Arguments:
//   indicies       Vector of grid indices
//   nsample        Vector for the number of grid cells in each direction
//
// Programmer:  Cody Raskin
// Creation:    July 1, 2015
//
// ****************************************************************************

inline
int
latticeIndex(const vector<int>& indicies,
             const vector<int>& nsample) {
    if (nsample.size() == 1) {
        return indicies[0];
    } else if (nsample.size() == 2) {
        return indicies[1]*nsample[0] + indicies[0];
    } else {
        return (indicies[2]*nsample[1]*nsample[0] +
                indicies[1]*nsample[0] +
                indicies[0]);
    }
}

// ****************************************************************************
// Method:  Helper Functions (mainly linear algebra)
//
// Purpose: Linear algebra
//
// Programmer:  Cody Raskin
// Creation:    July 1, 2015
//
// ****************************************************************************
template<int Dim>
inline
double
determinant(tensorStruct& T) {
    double result;
    if(Dim==3)  result = (T.xx*T.yy*T.zz + T.xy*T.yz*T.zx + T.xz*T.yx*T.zy) - (T.xx*T.yz*T.zy + T.xy*T.yx*T.zz + T.xz*T.yy*T.zx);
    else        result = T.xx*T.yy - T.xy*T.yx;
    return result;
}

template<int Dim>
inline
tensorStruct
square(tensorStruct& T) {
    tensorStruct S;
    if(Dim==3)
    {
        S.xx = T.xx*T.xx + T.xy*T.yx + T.xz*T.zx;
        S.xy = T.xx*T.xy + T.xy*T.yy + T.xz*T.zy;
        S.xz = T.xx*T.xz + T.xy*T.yz + T.xz*T.zz;
        S.yx = T.yx*T.xx + T.yy*T.yx + T.yz*T.zx;
        S.yy = T.yx*T.xy + T.yy*T.yy + T.yz*T.zy;
        S.yz = T.yx*T.xz + T.yy*T.yz + T.yz*T.zz;
        S.zx = T.zx*T.xx + T.zy*T.yx + T.zz*T.zx;
        S.zy = T.zx*T.xy + T.zy*T.yy + T.zz*T.zy;
        S.zz = T.zx*T.xz + T.zy*T.yz + T.zz*T.zz;
    }
    else
    {
        S.xx = T.xx*T.xx + T.xy*T.yx;
        S.xy = T.xx*T.xy + T.xy*T.yy;
        S.yx = T.yx*T.xx + T.yy*T.yx;
        S.yy = T.yx*T.xy + T.yy*T.yy;
    }
    return S;
}

template<int Dim>
inline
vector<double>
tensorVectorMult(tensorStruct& T,vector<double>& v) {
    vector<double> result(Dim);
    if(Dim==3)
    {
        result[0] = v[0]*T.xx + v[1]*T.xy + v[2]*T.xz;
        result[1] = v[0]*T.yx + v[1]*T.yy + v[2]*T.yz;
        result[2] = v[0]*T.zx + v[1]*T.zy + v[2]*T.zz;
    }
    else
    {
        result[0] = v[0]*T.xx + v[1]*T.xy;
        result[1] = v[0]*T.yx + v[1]*T.yy;
    }
    return result;
}

template<int Dim>
inline
tensorStruct
tensorScalarMult(tensorStruct& T, double x) {
    tensorStruct M;
    M.xx = T.xx * x;
    M.xy = T.xy * x;
    M.xz = T.xz * x;
    M.yx = T.yx * x;
    M.yy = T.yy * x;
    M.yz = T.yz * x;
    M.zx = T.zx * x;
    M.zy = T.zy * x;
    M.zz = T.zz * x;
    return M;
}

template<int Dim>
inline
tensorStruct
tensorSum(tensorStruct& T1, tensorStruct& T2) {
    tensorStruct M;
    M.xx = T1.xx + T2.xx;
    M.xy = T1.xy + T2.xy;
    M.xz = T1.xz + T2.xz;
    M.yx = T1.yx + T2.yx;
    M.yy = T1.yy + T2.yy;
    M.yz = T1.yz + T2.yz;
    M.zx = T1.zx + T2.zx;
    M.zy = T1.zy + T2.zy;
    M.zz = T1.zz + T2.zz;
    return M;
}

template<int Dim>
inline
tensorStruct
tensorInverse(tensorStruct& T) {
    tensorStruct M;
    double Tdet = determinant<Dim>(T);
    if(Dim==3)
    {
        M.xx = (T.yy*T.zz - T.yz*T.zy)/Tdet;
        M.xy = (T.xz*T.zy - T.xy*T.zz)/Tdet;
        M.xz = (T.xy*T.yz - T.xz*T.yy)/Tdet;
        M.yx = (T.yz*T.zx - T.yx*T.zz)/Tdet;
        M.yy = (T.xx*T.zz - T.xz*T.zx)/Tdet;
        M.yz = (T.xz*T.yx - T.xx*T.yz)/Tdet;
        M.zx = (T.yx*T.zy - T.yy*T.zx)/Tdet;
        M.zy = (T.xy*T.zx - T.xx*T.zy)/Tdet;
        M.zz = (T.xx*T.yy - T.xy*T.yx)/Tdet;
    }
    else
    {
        M.xx = T.yy/Tdet;
        M.xy = -1.0*T.xy/Tdet;
        M.yx = -1.0*T.yx/Tdet;
        M.yy = T.xx/Tdet;
    }
    return M;
}

template<int Dim>
inline
double
dotP(vector<double>& v1, vector<double>& v2) {
    double result;
    for(int k=0;k<Dim;++k)
        result += v1[k]*v2[k];
    return result;
}

template<int Dim>
inline
double
vectorMag(vector<double>& v) {
    double result=0;
    for(int d=0;d<Dim;++d)
        result += v[d]*v[d];
    result = sqrt(result);
    return result;
}

template<int Dim>
inline
vector<double>
vectorScalarMult(vector<double>& v, double s) {
    vector<double> result(Dim);
    for(int i=0;i<Dim;++i)
        result[i] = v[i]*s;
    return result;
}

template<int Dim>
inline
vector<double>
vectorSum(vector<double>& v1, vector<double>& v2) {
    vector<double> result(Dim);
    for(int i=0;i<Dim;++i)
        result[i] = v1[i] + v2[i];
    return result;
}

template<int Dim>
inline
tensorStruct
outerProduct(vector<double>& v1, vector<double>& v2) {
    tensorStruct M;
    M.xx = v1[0]*v2[0];
    M.xy = v1[0]*v2[1];
    M.yx = v1[1]*v2[0];
    M.yy = v1[1]*v2[1];
    if(Dim==3)
    {
        M.xz = v1[0]*v2[2];
        M.yz = v1[1]*v2[2];
        M.zx = v1[2]*v2[0];
        M.zy = v1[2]*v2[1];
        M.zz = v1[2]*v2[2];
    }
    return M;
}

template<int Dim>
tensorStruct *
avtSPHResampleFilter::CreateTensorStruct(vtkDataArray* const input_support, const int i)
{
    tensorStruct* Hi = new tensorStruct;
    
    if(input_support->GetNumberOfComponents() == 1)
    {
        if (Dim==3)
        {
            Hi->xx = 1.0/(*input_support->GetTuple(i));
            Hi->xy = 0.0;
            Hi->xz = 0.0;
            Hi->yx = 0.0;
            Hi->yy = 1.0/(*input_support->GetTuple(i));
            Hi->yz = 0.0;
            Hi->zx = 0.0;
            Hi->zy = 0.0;
            Hi->zz = 1.0/(*input_support->GetTuple(i));
        }
        else
        {
            Hi->xx = 1.0/(*input_support->GetTuple(i));
            Hi->xy = 0.0;
            Hi->yx = 0.0;
            Hi->yy = 1.0/(*input_support->GetTuple(i));
        }
    }
    else
    {
        if (Dim==3)
        {
            Hi->xx = (*input_support->GetTuple(i));
            Hi->xy = (*(input_support->GetTuple(i)+1));
            Hi->xz = (*(input_support->GetTuple(i)+2));
            Hi->yx = (*(input_support->GetTuple(i)+3));
            Hi->yy = (*(input_support->GetTuple(i)+4));
            Hi->yz = (*(input_support->GetTuple(i)+5));
            Hi->zx = (*(input_support->GetTuple(i)+6));
            Hi->zy = (*(input_support->GetTuple(i)+7));
            Hi->zz = (*(input_support->GetTuple(i)+8));
        }
        else
        {
            Hi->xx = (*input_support->GetTuple(i));
            Hi->xy = (*(input_support->GetTuple(i)+1));
            Hi->yx = (*(input_support->GetTuple(i)+3));
            Hi->yy = (*(input_support->GetTuple(i)+4));
        }
    }
    
    return Hi;
}

// ****************************************************************************
//  Method:      avtSPHResampleFilter::SampleNMS
//
//  Purpose:
//      Sample a given particle dataset to a rectilinear grid. This method
//      is called for serial execution only.
//
//  Arguments:
//      scalarValues    Vector of doubles holding the sampled data in a flat array
//
//  Programmer: Cody Raskin
//  Creation:   Fri Dec 5 13:50:31 PST 2014
//
//  Modifications:
//
//      Kevin Griffin, Thu Mar 24 18:28:09 PDT 2016
//      Added check for datasets = 0 and for NULL PointData to prevent Engine
//      crash.
//
//      Kevin Griffin, Wed May 25 14:38:01 PDT 2016
//      Cleaned-up the code, made variable names more meaningful, and removed
//      code designed to run in parallel.
//
// ****************************************************************************
template<int Dim>
void
avtSPHResampleFilter::SampleNMS(vector<double>& scalarValues)
{
    int timerHandle = visitTimer->StartTimer();
    std::string varname = resampleVarName;
    int nsets = 0;
    
    avtDataTree_p input_tree = GetInputDataTree();   /* get input data tree to obtain datasets */
    
    //    int num_input_vars = GetInput()->GetInfo().GetAttributes().GetNumberOfVariables();
    //    for most plots -- we expect at least 2 vars, one input var and one support tensor.
    vtkDataSet **data_sets = input_tree->GetAllLeaves(nsets);
    
    /// get the var we want to resample
    // make sure the name passed isn't simply our mesh name
    bool vars_to_resample = (GetInput()->GetInfo().GetAttributes().GetMeshname() != varname);
    
    // creating the lattice
    vector<double> latticeMin(Dim);
    vector<double> latticeMax(Dim);
    vector<int> gridCells(Dim); // X, Y, Z Grid Cells
    
    gridCells[0] = atts.GetXnum();
    gridCells[1] = atts.GetYnum();
    latticeMin[0] = atts.GetMinX();
    latticeMin[1] = atts.GetMinY();
    latticeMax[0] = atts.GetMaxX();
    latticeMax[1] = atts.GetMaxY();
    
    int totalCells = gridCells[0]*gridCells[1];
    
    if(Dim==3)
    {
        latticeMin[2] = atts.GetMinZ();
        latticeMax[2] = atts.GetMaxZ();
        gridCells[2] = atts.GetZnum();
        totalCells = totalCells * gridCells[2];
    }
    
    // Compute the step size.
    vector<double> stepXYZ = stepSize<Dim>(latticeMin, latticeMax, gridCells);
    
    debug5 << "Sampling to " << totalCells << " cells..." << endl;
    
    vector<int> partIdx = GetParticipatingIndices(data_sets, nsets);
    nsets = partIdx.size();
    
    if(vars_to_resample && (nsets > 0))
    {
        // number of particles
        vtkDataSet *dset = data_sets[partIdx[0]];
        vtkAppendPolyData *appender = NULL;
        
        if(nsets > 1)
        {
            appender = vtkAppendPolyData::New();
            
            for(int i=0; i<nsets; i++)
            {
                // Merge data sets
                appender->AddInputData(dynamic_cast<vtkPolyData *>(data_sets[partIdx[i]]));
            }
            
            appender->Update();
            dset = appender->GetOutput();
        }
        
        vtkIdType npart = dset->GetNumberOfPoints();;
        if (Dim==3)
        {
            debug5 << npart << " particles to " << gridCells[0] << "x" << gridCells[1] << "x" << gridCells[2] << " cells" << endl;
        }
        else
        {
            debug5 << npart << " particles to " << gridCells[0] << "x" << gridCells[1] << " cells" << endl;
        }
        
        
        
        // Prepare the output
        vector<double> m0(totalCells, 0.0);  // moments of data points
        vector<vector<double> > m1; // data points?
        vector<tensorStruct> m2(totalCells);
        vector<double> A(totalCells, 0.0);   // correction terms
        vector<vector<double> > B;  // correction terms
        
        scalarValues.resize(totalCells, 0.0);
        
        for(int k=0;k<totalCells;++k)
        {
            m1.push_back(vector<double>(Dim,0.0));
            B.push_back(vector<double>(Dim,0.0));
        }
        
        vtkDataArray *input_var = NULL;
        vtkDataArray *input_support = NULL;
        vtkDataArray *input_weight = NULL;
        
        input_var = dset->GetPointData()->GetArray(varname.c_str());
        
        if(input_var != NULL) {
            input_support = dset->GetPointData()->GetArray(supportVarName.c_str());
            input_weight  = dset->GetPointData()->GetArray(weightVarName.c_str());
        } else {
            input_var = dset->GetCellData()->GetArray(varname.c_str());
            input_support = dset->GetCellData()->GetArray(supportVarName.c_str());
            input_weight  = dset->GetCellData()->GetArray(weightVarName.c_str());
        }
        
        CheckInputVariables(input_var, input_support, input_weight);
        
        if(atts.GetRK())    // Don't need to compute the moments and correction values if RPK corrections is not selected
        {
            for(int i=0; i<npart; i++)
            {
                tensorStruct* Hi = CreateTensorStruct<Dim>(input_support, i);
                
                double Hdet = determinant<Dim>(*Hi);
                tensorStruct M = square<Dim>(*Hi);
                
                double weighti = *(input_weight->GetTuple(i));
                
                // Calculate smoothing length
                vector<double> extent(Dim);
                if(Dim==3)
                {
                    extent[0] = KERNEL_EXTENT/Hdet*sqrt(M.yy*M.zz - M.yz*M.zy);
                    extent[1] = KERNEL_EXTENT/Hdet*sqrt(M.xx*M.zz - M.xz*M.zx);
                    extent[2] = KERNEL_EXTENT/Hdet*sqrt(M.xx*M.yy - M.xy*M.yx);
                }
                else
                {
                    extent[0] = KERNEL_EXTENT/Hdet*sqrt(M.yy);
                    extent[1] = KERNEL_EXTENT/Hdet*sqrt(M.xx);
                }
                
                double ri[3] = {0, 0, 0};
                dset->GetPoint(i, &ri[0]);
                
                // Compute the set of lattice ids this node contributes to.
                const vector< vector<int> > ids = latticePoints<Dim>(ri, extent, latticeMin, stepXYZ, gridCells);
                
                // Iterate over the lattice ids.
                for (vector< vector<int> >::const_iterator itr = ids.begin(); itr != ids.end(); ++itr)
                {
                    vector<double> rj = latticePosition<Dim>(*itr, latticeMin, stepXYZ);    // Cell Center (xyz coords)
                    const int j = latticeIndex(*itr, gridCells);    // 1D Lattice Index
                    
                    // Contribution of the SPH point to this position.
                    vector<double> rij(Dim);
                    for(int k=0;k<Dim;++k)
                        rij[k] = ri[k]-rj[k];
                    
                    vector<double> etaVec(Dim);
                    etaVec              = tensorVectorMult<Dim>(*Hi,rij);
                    
                    const double etai   = vectorMag<Dim>(etaVec);
                    
                    //= (Hi*rij).magnitude();
                    const double Wi     = avtSPHResampleFilter::kernelValue<Dim>(etai, Hdet);
                    //                    const double gradWi = avtSPHResampleFilter::kernelGradValue<Dim>(etai, Hdet);
                    
                    // zeroth moment
                    const double wwi    = Wi*weighti;
                    m0[j]               += wwi;
                    //gradm0[j]           += gradWi*weighti;
                    
                    // first moment
                    vector<double> m1i(Dim);
                    m1i                 = vectorScalarMult<Dim>(rij,wwi);
                    m1[j]               = vectorSum<Dim>(m1[j],m1i);
                    /*
                     vector<double> vgradWi(Dim);
                     for(int k=0;k<Dim;++k) vgradWi[k] = gradWi;
                     vector<double> rji  = vectorScalarMult<Dim>(rij,-1.0);
                     tensorStruct gm1i   = outerProduct<Dim>(rji,vgradWi);
                     gm1i.xx             += Wi;
                     gm1i.yy             += Wi;
                     gm1i.zz             += Wi;
                     gm1i                = tensorScalarMult<Dim>(gm1i,weighti);
                     gradm1[j]           = tensorSum<Dim>(gradm1[j],gm1i);
                     */
                    
                    //second moment
                    tensorStruct thpt;
                    thpt                = outerProduct<Dim>(rij,rij);
                    tensorStruct m2i;
                    m2i                 = tensorScalarMult<Dim>(thpt,wwi);
                    m2[j]               = tensorSum<Dim>(m2[j],m2i);
                    
                    /*
                     printf("\n\n< ");
                     for(int k=0;k<Dim;++k)
                     printf("%f ",rj[k]);
                     printf(">\n");
                     
                     printf("i=%d\n",i);
                     printf("rij=<%3.3f,%3.3f>\n",rij[0],rij[1]);
                     printf("etai=%3.2f Wi=%3.2f wwi=%3.2f\n",etai,Wi,wwi);
                     
                     printf("m1=<%3.3f,%3.3f>\n",m1i[0],m1i[1]);
                     
                     printf("m2=[[%3.3f,%3.3f][%3.3f,%3.3f]]\n",m2i.xx,m2i.xy,m2i.yx,m2i.yy);
                     */
                    
                }
                
                delete Hi;
            }
        
            // compute A & B on the lattice points
            for(int i=0; i<totalCells;++i)
            {
                double m2det            = determinant<Dim>(m2[i]);
                //double m2det = 0;
                tensorStruct m2inv;
                if(fabs(m2det) < 1.0e-10) m2det = 0;
                if(fabs(m2det) > 0)
                m2inv               = tensorInverse<Dim>(m2[i]);
            
                //printf("m2inv=[[%3.2f %3.2f][%3.2f %3.2f]]\n",m2inv.xx,m2inv.xy,m2inv.yx,m2inv.yy);
            
                vector<double> m2invm1(Dim);
                m2invm1                 = tensorVectorMult<Dim>(m2inv,m1[i]);
            
                //printf("m2invm1=[%3.3f,%3.3f]\n",m2invm1[0],m2invm1[1]);
                //double Ainv             = m0[i] - dotP<Dim>(m2invm1,m1[i]);
                double Ainv             = m0[i];
                A[i]                    = (fabs(Ainv)>1e-10 ? 1.0/Ainv : 1.0);
                B[i]                    = vectorScalarMult<Dim>(m2invm1,-1.0);
            }
        } // End if(atts.GetRK())
        
        for(int i=0; i<npart;++i)
        {
            tensorStruct* Hi = CreateTensorStruct<Dim>(input_support, i);

            double Hdet = determinant<Dim>(*Hi);
            tensorStruct M = square<Dim>(*Hi);
            
            double weighti = *input_weight->GetTuple(i);
            
            vector<double> extent(Dim);
            if(Dim==3)
            {
                extent[0] = KERNEL_EXTENT/Hdet*sqrt(M.yy*M.zz - M.yz*M.zy);
                extent[1] = KERNEL_EXTENT/Hdet*sqrt(M.xx*M.zz - M.xz*M.zx);
                extent[2] = KERNEL_EXTENT/Hdet*sqrt(M.xx*M.yy - M.xy*M.yx);
            }
            else
            {
                extent[0] = KERNEL_EXTENT/Hdet*sqrt(M.yy);
                extent[1] = KERNEL_EXTENT/Hdet*sqrt(M.xx);
            }
            
            double ri[3] = {0, 0, 0};
            dset->GetPoint(i, &ri[0]);
            
            // Compute the set of lattice ids this node contributes to.
            const vector< vector<int> > ids = latticePoints<Dim>(ri, extent, latticeMin, stepXYZ, gridCells);
            
            double field = *input_var->GetTuple(i);
            
            // Iterate over the lattice points.
            for (vector< vector<int> >::const_iterator itr = ids.begin();
                 itr != ids.end();
                 ++itr)
            {
                // Compute the sample position and index into the lattice sample array.
                vector<double> rj = latticePosition<Dim>(*itr, latticeMin, stepXYZ);
                const int j = latticeIndex(*itr, gridCells);
                
                // Contribution of the SPH point to this position.
                vector<double> rij(Dim);
                for(int k=0;k<Dim;++k)
                    rij[k] = ri[k]-rj[k];
                vector<double> etaVec(Dim);
                etaVec            = tensorVectorMult<Dim>(*Hi,rij);
                
                const double etai = vectorMag<Dim>(etaVec);
                
                //= (Hi*rij).magnitude();
                const double Wi = avtSPHResampleFilter::kernelValue<Dim>(etai, Hdet);
                const double thpt = weighti*Wi;
                //scalarValues[j] += thpt * field;
                if(atts.GetRK())
                    //scalarValues[j] += thpt * field * A[j] * (1.0+dotP<Dim>(B[j],rij));
                    scalarValues[j] += thpt * field * A[j];
                else
                    scalarValues[j] += thpt * field;
                //printf("SPH=%3.2e RK=%3.2e\n",thpt * field,thpt * field * (1.0 + A[j] * (1.0+dotP<Dim>(B[j],rij))));
            }
            
            delete Hi;
        }
        
        if(appender != NULL)
        {
            appender->Delete();
        }
    }
    else
    {
        scalarValues.resize(totalCells, 0.0);
    }
    
    visitTimer->StopTimer(timerHandle, "avtSPHResample::SampleNMS");
}

// ****************************************************************************
//  Method:      avtSPHResampleFilter::Sample
//
//  Purpose:
//      Sample a given particle dataset to a rectilinear grid.
//
//  Arguments:
//      data    Vector of of SPH data. Currently this vector will contain at
//              most one item.

//
//  Programmer: Cody Raskin
//  Creation:   Fri Dec 5 13:50:31 PST 2014
//
//  Modifications:
//
//      Kevin Griffin, Thu Mar 24 18:28:09 PDT 2016
//      Added check for datasets = 0 and for NULL PointData to prevent Engine
//      crash.
//
//      Kevin Griffin, Tue Jun 21 13:09:21 PDT 2016
//      Reimplemented to work better in parallel and scale in memory.
//
//      Kevin Griffin, Wed Aug 17 16:17:16 PDT 2016
//      Extended the dataset bounds, if needed, to match the user supplied
//      min/max bounds. Fixed the issue with holes appearing in the final
//      image.
//
// ****************************************************************************
template<int Dim>
void
avtSPHResampleFilter::Sample(vector<sphData>& data)
{
    int timerHandle = visitTimer->StartTimer();
    
    int nsets = 0;
    vector<int> partIdx;
    avtDataTree_p input_tree = GetInputDataTree();
    
    vtkDataSet **data_sets = input_tree->GetAllLeaves(nsets);
    
    // Global Lattice Info
    vector<double> globalLatticeMin(Dim);
    vector<double> globalLatticeMax(Dim);
    vector<int> globalGridCells(Dim);
    
    globalGridCells[0] = atts.GetXnum();
    globalGridCells[1] = atts.GetYnum();
    globalLatticeMin[0] = atts.GetMinX();
    globalLatticeMin[1] = atts.GetMinY();
    globalLatticeMax[0] = atts.GetMaxX();
    globalLatticeMax[1] = atts.GetMaxY();
    
    if(Dim==3)
    {
        globalLatticeMin[2] = atts.GetMinZ();
        globalLatticeMax[2] = atts.GetMaxZ();
        globalGridCells[2] = atts.GetZnum();
    }
    
    // Local Lattice Info
    vector<double> latticeMin(Dim);
    vector<double> latticeMax(Dim);
    vector<int> gridCells(Dim);
    
    int totalCells  = 0;
    vector<double> stepXYZ = GetStepSizeXYZ();
    
    // Determine participation and merge datasets
    if(nsets > 0)
    {
        partIdx = GetParticipatingIndices(data_sets, nsets);
        nsets = partIdx.size();
    }
    
    bool isParticipating = nsets > 0;
    vtkDataSet *dset = isParticipating ? data_sets[partIdx[0]] : NULL;
    vtkAppendPolyData *appender = NULL;
    
    if(nsets > 1)
    {
        appender = vtkAppendPolyData::New();
        
        for(int i=0; i<nsets; i++)
        {
            // Merge data sets
            appender->AddInputData(dynamic_cast<vtkPolyData *>(data_sets[partIdx[i]]));
        }
        
        appender->Update();
        dset = appender->GetOutput();
    }
    
    if(atts.GetRK())
    {
        if(isParticipating)
        {
            for(int i=0; i<1; i++)  // For now I'm merging the datasets into 1, change if decided to process each domain individually
            {
                // Number of particles
                vtkIdType npart = dset->GetNumberOfPoints();
                
                double bounds[6];
                dset->GetBounds(bounds);
                
                ExtendBoundsIfNeeded(bounds);
                
                debug5 << "Bounds " << bounds[0] << ", " << bounds[1] << ", " << bounds[2] << ", " << bounds[3] << ", " << bounds[4] << ", " << bounds[5] << endl;
                
                GetDimensions(gridCells, latticeMin, latticeMax, stepXYZ, bounds);
                
                if(Dim == 2)
                {
                    totalCells = gridCells[0]*gridCells[1];
                    debug5 << npart << " particles to " << gridCells[0] << "x" << gridCells[1] << " cells" << endl;
                }
                else
                {
                    totalCells = gridCells[0]*gridCells[1]*gridCells[2];
                    debug5 << npart << " particles to " << gridCells[0] << "x" << gridCells[1] << "x" << gridCells[2] << " cells" << endl;
                }
                
                debug5 << "Sampling to " << totalCells << " cells..." << endl;
                
                sphData currentData = sphData(totalCells, Dim);
                currentData.dset = dset;
                currentData.bounds[0] = bounds[0];
                currentData.bounds[1] = bounds[1];
                currentData.bounds[2] = bounds[2];
                currentData.bounds[3] = bounds[3];
                currentData.bounds[4] = bounds[4];
                currentData.bounds[5] = bounds[5];
                
                currentData.indices[0] = max(0, int(floor((latticeMin[0] - globalLatticeMin[0]) / stepXYZ[0])));
                currentData.indices[1] = currentData.indices[0] + gridCells[0];
                
                currentData.indices[2] = max(0, int(floor((latticeMin[1] - globalLatticeMin[1]) / stepXYZ[1])));;
                currentData.indices[3] = currentData.indices[2] + gridCells[1];
                
                if(Dim == 3)
                {
                    currentData.indices[4] = max(0, int(floor((latticeMin[2] - globalLatticeMin[2]) / stepXYZ[2])));;
                    currentData.indices[5] = currentData.indices[4] + gridCells[2];
                    
                    debug5 << "Indices " << currentData.indices[0] << ", " << currentData.indices[1] << ", " << currentData.indices[2] << ", " << currentData.indices[3] << ", " << currentData.indices[4] << ", " << currentData.indices[5] << endl;
                }
                else
                {
                    debug5 << "Indices " << currentData.indices[0] << ", " << currentData.indices[1] << ", " << currentData.indices[2] << ", " << currentData.indices[3] << endl;
                }
                
                // Get Input Variables
                vtkDataArray *input_var = dset->GetPointData()->GetArray(resampleVarName.c_str());
                vtkDataArray *input_support = NULL;
                vtkDataArray *input_weight = NULL;
                
                if(input_var != NULL) {
                    input_support = dset->GetPointData()->GetArray(supportVarName.c_str());
                    input_weight  = dset->GetPointData()->GetArray(weightVarName.c_str());
                } else {
                    input_var = dset->GetCellData()->GetArray(resampleVarName.c_str());
                    input_support = dset->GetCellData()->GetArray(supportVarName.c_str());
                    input_weight  = dset->GetCellData()->GetArray(weightVarName.c_str());
                }
                
                CheckInputVariables(input_var, input_support, input_weight);
                
                // Calculate moments
                for(int ii=0; ii<npart; ii++)
                {
                    tensorStruct *Hi = CreateTensorStruct<Dim>(input_support, ii);
                    
                    double Hdet = determinant<Dim>(*Hi);
                    tensorStruct M = square<Dim>(*Hi);
                    double weighti = *(input_weight->GetTuple(ii));
                    
                    // Calculate Smoothing Length
                    vector<double> extent(Dim);
                    extent[0] = KERNEL_EXTENT/Hdet*sqrt(M.yy);
                    extent[1] = KERNEL_EXTENT/Hdet*sqrt(M.xx);
                    
                    if(Dim == 3)
                    {
                        extent[0] = KERNEL_EXTENT/Hdet*sqrt(M.yy*M.zz - M.yz*M.zy);
                        extent[1] = KERNEL_EXTENT/Hdet*sqrt(M.xx*M.zz - M.xz*M.zx);
                        extent[2] = KERNEL_EXTENT/Hdet*sqrt(M.xx*M.yy - M.xy*M.yx);
                    }
                    
                    double ri[3] = {0,0,0};
                    dset->GetPoint(ii, &ri[0]);
                    
                    // Compute the set of lattice ids this node contributes to.
                    const vector< vector<int> > ids = latticePoints<Dim>(ri, extent, globalLatticeMin, stepXYZ, globalGridCells);
                    ghostCell temp(Dim);
                    
                    // Iterate over the lattice ids.
                    for (vector< vector<int> >::const_iterator itr = ids.begin(); itr != ids.end(); ++itr)
                    {
                        vector<double> rj = latticePosition<Dim>(*itr, globalLatticeMin, stepXYZ);    // Cell Center (xyz coords)
                        const int j = latticeIndex(*itr, globalGridCells);    // 1D Lattice Index
                        int localj = GetLocalIndex(currentData, j);
                        temp.globalLatticeIdx = j;
                        
                        // Contribution of the SPH point to this position.
                        vector<double> rij(Dim);
                        for(int k=0;k<Dim;++k)
                        {
                            rij[k] = ri[k]-rj[k];
                        }
                        
                        vector<double> etaVec(Dim);
                        etaVec = tensorVectorMult<Dim>(*Hi,rij);
                        
                        const double etai = vectorMag<Dim>(etaVec);
                        const double Wi = avtSPHResampleFilter::kernelValue<Dim>(etai, Hdet);
                        
                        // zeroth moment
                        const double wwi = Wi * weighti;
                        temp.m0 = wwi;
                        
                        // first moment
                        temp.m1 = vectorScalarMult<Dim>(rij,wwi);
                        
                        // second moment
                        tensorStruct thpt;
                        thpt = outerProduct<Dim>(rij,rij);
                        temp.m2 = tensorScalarMult<Dim>(thpt,wwi);
                        
                        if(localj >= 0)
                        {
                            currentData.globalLatticeIndices[localj] = temp.globalLatticeIdx;
                            
                            currentData.m0[localj] += temp.m0;
                            currentData.m1[localj] = vectorSum<Dim>(currentData.m1[localj], temp.m1);
                            currentData.m2[localj] = tensorSum<Dim>(currentData.m2[localj], temp.m2);
                        }
                        else
                        {
                            int idx = GetGhostCellIndex(currentData.ghostCells, j);
                            
                            if(idx >= 0)
                            {
                                currentData.ghostCells[idx].m0 += temp.m0;
                                currentData.ghostCells[idx].m1 = vectorSum<Dim>(currentData.ghostCells[idx].m1, temp.m1);
                                currentData.ghostCells[idx].m2 = tensorSum<Dim>(currentData.ghostCells[idx].m2, temp.m2);
                            }
                            else
                            {
                                ghostCell gc(Dim);
                                
                                gc.globalLatticeIdx = temp.globalLatticeIdx;
                                gc.m0 = temp.m0;
                                gc.m1 = temp.m1;
                                gc.m2 = temp.m2;
                                
                                currentData.ghostCells.push_back(gc);
                            }
                        }
                    }
                    
                    delete Hi;
                } // for(npart)
                debug5 << "Total Cells: " << currentData.totalCells << " Ghost Cells: " << currentData.ghostCells.size() << endl;
                data.push_back(currentData);
            }   // End For(nsets)
        }   // isParticipating
        else
        {
            double bounds[6] = {0,0,0,0,0,0};
            ExtendBoundsIfNeeded(bounds);
        }
        
         visitTimer->StopTimer(timerHandle, "avtSPHResample::Sample - stop 1");
        
        // Synchronize moments
        if(data.size() > 0)
        {
            SynchMoments(data[0]);
        }
        else
        {
            sphData fakeData(0, Dim);
            SynchMoments(fakeData);
        }
        
        timerHandle = visitTimer->StartTimer();
        
        // Compute A and B Correction Values
        for(int i=0; i<data.size(); i++)
        {
            for(int j=0; j<data[i].totalCells; j++)
            {
                double m2det            = determinant<Dim>(data[i].m2[j]);
                //double m2det = 0;
                tensorStruct m2inv;
                if(fabs(m2det) < 1.0e-10)
                {
                    m2det = 0;
                }
                
                if(fabs(m2det) > 0)
                {
                    m2inv               = tensorInverse<Dim>(data[i].m2[j]);
                }
                
                vector<double> m2invm1(Dim);
                m2invm1                 = tensorVectorMult<Dim>(m2inv,data[i].m1[j]);
                double Ainv             = data[i].m0[j];
                data[i].A[j]                    = (fabs(Ainv)>1e-10 ? 1.0/Ainv : 1.0);
                data[i].B[j]                    = vectorScalarMult<Dim>(m2invm1,-1.0);
            }
            
            for(int j=0; j<data[i].ghostCells.size(); j++)
            {
                ghostCell *gc = &(data[i].ghostCells[j]);
                double m2det            = determinant<Dim>(gc->m2);
                //double m2det = 0;
                tensorStruct m2inv;
                if(fabs(m2det) < 1.0e-10)
                {
                    m2det = 0;
                }
                
                if(fabs(m2det) > 0)
                {
                    m2inv               = tensorInverse<Dim>(gc->m2);
                }
                
                vector<double> m2invm1(Dim);
                m2invm1                 = tensorVectorMult<Dim>(m2inv,gc->m1);
                double Ainv             = gc->m0;
                gc->A                    = (fabs(Ainv)>1e-10 ? 1.0/Ainv : 1.0);
                gc->B                    = vectorScalarMult<Dim>(m2invm1,-1.0);
            }
        }
        
        // Compute Scalar Values
        for(int i=0; i<data.size(); i++)
        {
            vtkDataSet *dset = data[i].dset;
            vtkIdType npart = dset->GetNumberOfPoints();
//            double *bounds = data[i].bounds;
            
            // Get Input Variables
            vtkDataArray *input_var = dset->GetPointData()->GetArray(resampleVarName.c_str());
            vtkDataArray *input_support = NULL;
            vtkDataArray *input_weight = NULL;
            
            if(input_var != NULL) {
                input_support = dset->GetPointData()->GetArray(supportVarName.c_str());
                input_weight  = dset->GetPointData()->GetArray(weightVarName.c_str());
            } else {
                input_var = dset->GetCellData()->GetArray(resampleVarName.c_str());
                input_support = dset->GetCellData()->GetArray(supportVarName.c_str());
                input_weight  = dset->GetCellData()->GetArray(weightVarName.c_str());
            }
            
            for(int ii=0; ii<npart; ii++)
            {
                tensorStruct* Hi = CreateTensorStruct<Dim>(input_support, ii);
                
                double Hdet = determinant<Dim>(*Hi);
                tensorStruct M = square<Dim>(*Hi);
                
                double weighti = *input_weight->GetTuple(ii);
                
                vector<double> extent(Dim);
                if(Dim==3)
                {
                    extent[0] = KERNEL_EXTENT/Hdet*sqrt(M.yy*M.zz - M.yz*M.zy);
                    extent[1] = KERNEL_EXTENT/Hdet*sqrt(M.xx*M.zz - M.xz*M.zx);
                    extent[2] = KERNEL_EXTENT/Hdet*sqrt(M.xx*M.yy - M.xy*M.yx);
                }
                else
                {
                    extent[0] = KERNEL_EXTENT/Hdet*sqrt(M.yy);
                    extent[1] = KERNEL_EXTENT/Hdet*sqrt(M.xx);
                }
                
                double ri[3] = {0, 0, 0};
                dset->GetPoint(ii, &ri[0]);
                
                // Compute the set of lattice ids this node contributes to.
                const vector< vector<int> > ids = latticePoints<Dim>(ri, extent, globalLatticeMin, stepXYZ, globalGridCells);
                
                double field = *input_var->GetTuple(ii);
                
                // Iterate over the lattice points.
                for (vector< vector<int> >::const_iterator itr = ids.begin(); itr != ids.end(); ++itr)
                {
                    // Compute the sample position and index into the lattice sample array.
                    vector<double> rj = latticePosition<Dim>(*itr, globalLatticeMin, stepXYZ);
                    const int j = latticeIndex(*itr, globalGridCells);
                    int localj = GetLocalIndex(data[i], j);
                    
                    // Contribution of the SPH point to this position.
                    vector<double> rij(Dim);
                    for(int k=0;k<Dim;++k)
                    {
                        rij[k] = ri[k]-rj[k];
                    }
                    
                    vector<double> etaVec(Dim);
                    etaVec            = tensorVectorMult<Dim>(*Hi,rij);
                    
                    const double etai = vectorMag<Dim>(etaVec);
                    
                    //= (Hi*rij).magnitude();
                    const double Wi = avtSPHResampleFilter::kernelValue<Dim>(etai, Hdet);
                    const double thpt = weighti*Wi;
                    
                    if(localj >= 0)
                    {
                        data[i].scalarValues[localj] += thpt * field * data[i].A[localj];
                    }
                    else
                    {
                        int idx = GetGhostCellIndex(data[i].ghostCells, j);
                        
                        if(idx >= 0)
                        {
                            data[i].ghostCells[idx].value += thpt * field * data[i].ghostCells[idx].A;
                        }
                        else    // Should never get here
                        {
                            EXCEPTION1(VisItException, "Couldn't find a Ghost Cell that should exist");
                        }
                    }
                }
                delete Hi;
            }   // for(npart)
        } // for(data.size())
    } // End if(atts.GetRPK())
    else
    {
        if(isParticipating)
        {
            // Compute Scalar Values without RPK Corrections
            double bounds[6];
            dset->GetBounds(bounds);
            ExtendBoundsIfNeeded(bounds);
            debug5 << "Bounds " << bounds[0] << ", " << bounds[1] << ", " << bounds[2] << ", " << bounds[3] << ", " << bounds[4] << ", " << bounds[5] << endl;
            
            GetDimensions(gridCells, latticeMin, latticeMax, stepXYZ, bounds);
            totalCells = Dim == 3 ? gridCells[0]*gridCells[1]*gridCells[2] : gridCells[0]*gridCells[1];
            
            debug5 << "Sampling to " << totalCells << " cells..." << endl;
            
            sphData currentData = sphData(totalCells, Dim);
            currentData.dset = dset;
            currentData.bounds[0] = bounds[0];
            currentData.bounds[1] = bounds[1];
            currentData.bounds[2] = bounds[2];
            currentData.bounds[3] = bounds[3];
            currentData.bounds[4] = bounds[4];
            currentData.bounds[5] = bounds[5];
            
            currentData.indices[0] = max(0, int(floor((latticeMin[0] - globalLatticeMin[0]) / stepXYZ[0])));
            currentData.indices[1] = currentData.indices[0] + gridCells[0];
            
            currentData.indices[2] = max(0, int(floor((latticeMin[1] - globalLatticeMin[1]) / stepXYZ[1])));;
            currentData.indices[3] = currentData.indices[2] + gridCells[1];
            
            if(Dim == 3)
            {
                currentData.indices[4] = max(0, int(floor((latticeMin[2] - globalLatticeMin[2]) / stepXYZ[2])));;
                currentData.indices[5] = currentData.indices[4] + gridCells[2];
            }
            
            // Get Input Variables
            vtkDataArray *input_var = dset->GetPointData()->GetArray(resampleVarName.c_str());
            vtkDataArray *input_support = NULL;
            vtkDataArray *input_weight = NULL;
            
            if(input_var != NULL) {
                input_support = dset->GetPointData()->GetArray(supportVarName.c_str());
                input_weight  = dset->GetPointData()->GetArray(weightVarName.c_str());
            } else {
                input_var = dset->GetCellData()->GetArray(resampleVarName.c_str());
                input_support = dset->GetCellData()->GetArray(supportVarName.c_str());
                input_weight  = dset->GetCellData()->GetArray(weightVarName.c_str());
            }
            
            CheckInputVariables(input_var, input_support, input_weight);
            
            vtkIdType npart = dset->GetNumberOfPoints();
            
            // Compute Scalar Values
            for(int ii=0; ii<npart; ii++)
            {
                tensorStruct* Hi = CreateTensorStruct<Dim>(input_support, ii);
                
                double Hdet = determinant<Dim>(*Hi);
                tensorStruct M = square<Dim>(*Hi);
                
                double weighti = *input_weight->GetTuple(ii);
                
                vector<double> extent(Dim);
                if(Dim==3)
                {
                    extent[0] = KERNEL_EXTENT/Hdet*sqrt(M.yy*M.zz - M.yz*M.zy);
                    extent[1] = KERNEL_EXTENT/Hdet*sqrt(M.xx*M.zz - M.xz*M.zx);
                    extent[2] = KERNEL_EXTENT/Hdet*sqrt(M.xx*M.yy - M.xy*M.yx);
                }
                else
                {
                    extent[0] = KERNEL_EXTENT/Hdet*sqrt(M.yy);
                    extent[1] = KERNEL_EXTENT/Hdet*sqrt(M.xx);
                }
                
                double ri[3] = {0, 0, 0};
                dset->GetPoint(ii, &ri[0]);
                
                // Compute the set of lattice ids this node contributes to.
                const vector< vector<int> > ids = latticePoints<Dim>(ri, extent, globalLatticeMin, stepXYZ, globalGridCells);
                
                double field = *input_var->GetTuple(ii);
                
                // Iterate over the lattice points.
                for (vector< vector<int> >::const_iterator itr = ids.begin();
                     itr != ids.end();
                     ++itr)
                {
                    // Compute the sample position and index into the lattice sample array.
                    vector<double> rj = latticePosition<Dim>(*itr, globalLatticeMin, stepXYZ);
                    const int j = latticeIndex(*itr, globalGridCells);
                    int localj = GetLocalIndex(currentData, j);
                    
                    // Contribution of the SPH point to this position.
                    vector<double> rij(Dim);
                    for(int k=0;k<Dim;++k)
                    {
                        rij[k] = ri[k]-rj[k];
                    }
                    
                    vector<double> etaVec(Dim);
                    etaVec            = tensorVectorMult<Dim>(*Hi,rij);
                    
                    const double etai = vectorMag<Dim>(etaVec);
                    
                    //= (Hi*rij).magnitude();
                    const double Wi = avtSPHResampleFilter::kernelValue<Dim>(etai, Hdet);
                    const double thpt = weighti*Wi;
                    
                    if(localj >= 0)
                    {
                        currentData.globalLatticeIndices[localj] = j;
                        currentData.scalarValues[localj] += thpt * field;
                    }
                    else
                    {
                        int idx = GetGhostCellIndex(currentData.ghostCells, j);
                        
                        if(idx >= 0)
                        {
                            currentData.ghostCells[idx].value += thpt * field;
                        }
                        else
                        {
                            ghostCell gc(Dim);
                            
                            gc.globalLatticeIdx = j;
                            gc.value += thpt * field;
                            
                            currentData.ghostCells.push_back(gc);
                        }
                    }
                }
                
                delete Hi;
            }   // for(npart)
            data.push_back(currentData);
        }   // End isParticipating
        else
        {
            double bounds[6] = {0,0,0,0,0,0};
            ExtendBoundsIfNeeded(bounds);
        }
    } // End Else
    
    visitTimer->StopTimer(timerHandle, "avtSPHResample::Sample - stop 2");
    
    // Synchronize scalar values
    if(data.size() > 0)
    {
        SynchScalars(data[0]);
    }
    else
    {
        sphData fakeData(0, Dim);
        SynchScalars(fakeData);
    }
    
    if(appender != NULL)
    {
        appender->Delete();
    }
    
//    visitTimer->StopTimer(timerHandle, "avtSPHResample::Sample");
}

// ****************************************************************************
//  Method: avtSPHResampleFilter::Execute
//
//  Purpose:
//      Sends the specified input and output through the SPHResample filter.
//
//  Returns:       The output dataset.
//
//  Programmer: harrison37 -- generated by xml2avt
//  Creation:   Fri Dec 5 13:50:31 PST 2014
//
//  Modifications:
//  Kevin Griffin, Mon Jan 4 17:51:20 PST 2016
//  Added a check for cell-centered data.
//
//  Kevin Griffin, Tue Jun 21 13:09:21 PDT 2016
//  Reimplemented to work better in parallel and scale in memory
//
// ****************************************************************************
void
avtSPHResampleFilter::Execute()
{
#ifndef PARALLEL
    ExecuteNMS();
    return;
#endif
    
    vector<sphData> vectorOfSPHData;
    
    avtDataAttributes &inAtts = GetInput()->GetInfo().GetAttributes();
    nDim = inAtts.GetSpatialDimension();
        
    if(nDim==3)
    {
        avtSPHResampleFilter::Sample<3>(vectorOfSPHData);
    }
    else
    {
        avtSPHResampleFilter::Sample<2>(vectorOfSPHData);
    }

    int size = vectorOfSPHData.size();
    if(size > 0)
    {
        int sampledPointsSize = vectorOfSPHData[0].scalarValues.size();
        
        // create output data array
        vtkDoubleArray *out_var = vtkDoubleArray::New();
        out_var->SetName(resampleVarName.c_str());
        out_var->SetNumberOfTuples(sampledPointsSize);
        out_var->SetNumberOfComponents(1);
        
        for(int j=0; j<sampledPointsSize; j++)
        {
            out_var->SetTuple1(j, vectorOfSPHData[0].scalarValues[j]);
        }
        
        vector<double> steps = GetStepSizeXYZ();
        vtkDataSet *out_dset = CreateOutputGrid(vectorOfSPHData[0].bounds, steps);
        
        out_dset->GetCellData()->AddArray(out_var);
        out_dset->GetCellData()->SetActiveScalars(out_var->GetName());
        out_var->Delete();
        
        avtDataTree_p out_tree = new avtDataTree(out_dset, 0);
        SetOutputDataTree(out_tree);
        out_dset->Delete();
    }
    else
    {
        debug5 << "Sending dummy data" << endl;
        avtDataTree_p dummy = new avtDataTree();
        SetOutputDataTree(dummy);
    }
}

// ****************************************************************************
//  Method: avtSPHResampleFilter::ExecuteNMS
//
//  Purpose:
//      Sends the specified input and output through the SPHResample filter.
//
//
//  Returns:       The output dataset.
//
//  Programmer: Kevin Griffin
//  Creation:   Fri Dec 5 13:50:31 PST 2014
//
//  Modifications:
//  Kevin Griffin, Mon Jan 4 17:51:20 PST 2016
//  Added a check for cell-centered data.
//
//  Kevin Griffin, Wed May 25 14:38:01 PDT 2016
//  Removed the check for cell-centered data. Cleaned-up the code.
//
// ****************************************************************************
void
avtSPHResampleFilter::ExecuteNMS()
{
    vector<double> vectorofScalars;
    
    avtDataAttributes &inAtts = GetInput()->GetInfo().GetAttributes();
    nDim = inAtts.GetSpatialDimension();
    
    if(nDim==3)
    {
        avtSPHResampleFilter::SampleNMS<3>(vectorofScalars);
    }
    else
    {
        avtSPHResampleFilter::SampleNMS<2>(vectorofScalars);
    }
    
    int sampledPointsSize = vectorofScalars.size();
    std::string varname = resampleVarName;
    
    // create output data array
    vtkDoubleArray *out_var = vtkDoubleArray::New();
    out_var->SetName(varname.c_str());
    out_var->SetNumberOfTuples(sampledPointsSize);
    out_var->SetNumberOfComponents(1);
        
    for(int i=0;i<sampledPointsSize;i++)
    {
        out_var->SetTuple1(i,vectorofScalars[i]);
    }
    
    double bounds[6] = {atts.GetMinX(), atts.GetMaxX(), atts.GetMinY(), atts.GetMaxY(), atts.GetMinZ(), atts.GetMaxZ()};
    
    vtkRectilinearGrid *out_dset = CreateOutputGrid(bounds, GetStepSizeXYZ());
    out_dset->GetCellData()->AddArray(out_var);
    out_dset->GetCellData()->SetActiveScalars(out_var->GetName());
    out_var->Delete();
    
    avtDataTree_p out_tree = new avtDataTree(out_dset, 0);
    out_dset->Delete();
    SetOutputDataTree(out_tree);
}

// ****************************************************************************
// Method:  avtSPHResampleFilter::GetParticipatingIndices
//
// Purpose: Determines the indices of the datasets that need to be processed.
//          Potentially, not all of the datasets will need to be processed if
//          the user selects a subset of the data to resample.
//
// Arguments:
//      dsets   The datasets that this rank has
//      nsets   The number of datasets
//
// Programmer:  Kevin Griffin
// Creation:    Tue Jun 21 16:37:00 PDT 2016
//
//  Modifications:
//      Kevin Griffin, Thu Dec 22 13:08:30 PST 2016
//      Commented out the filtering code. The filtering needs to use the
//      smoothing (influence) length in addtion to the bounds. The current
//      implementation caused too many particles to be removed which resulted in
//      the appearance of "holes" in the final image.
//
// ****************************************************************************
vector<int>
avtSPHResampleFilter::GetParticipatingIndices(vtkDataSet **dsets, const int nsets)
{
    vector<int> partIdx;
    /*bool good;
    double bounds[6];
    
    double latticeMinXYZ[3] = {atts.GetMinX(), atts.GetMinY(), nDim == 3 ? atts.GetMinZ() : 0};
    double latticeMaxXYZ[3] = {atts.GetMaxX(), atts.GetMaxY(), nDim == 3 ? atts.GetMaxZ() : 1};
    
    for(int j=0; j<nsets; j++)
    {
        good = true;
        dsets[j]->GetBounds(bounds);
        
        for(int i=0; i<nDim; i++)
        {
            int offset = i*2;
            if( (latticeMinXYZ[i] > bounds[offset+1]) || (latticeMaxXYZ[i] < bounds[offset]) )
            {
                good = false;
                break;
            }
        }
        
        if(good)
        {
            partIdx.push_back(j);
        }
    }*/
    for(int j=0; j<nsets; j++)
    {
        partIdx.push_back(j);
    }
    
    
    return partIdx;
}

// ****************************************************************************
// Method:  avtSPHResampleFilter::SynchScalars
//
// Purpose: Perform a sum reduction of the scalar values across all ranks.
//
// Arguments:
//      data  The SPH data containing this rank's scalar values
//
// Programmer:  Kevin Griffin
// Creation:    Thu Jun 16 19:07:57 PDT 2016
//
//  Modifications:
//
// ****************************************************************************
void
avtSPHResampleFilter::SynchScalars(sphData &data)
{
    int myRank = PAR_Rank();
    
    // Collect all the ghost cell data that need to be synchronized
    // Need to be done before values are updated
    vector<double> o_scalars;
    vector<int> o_latticeIndices;
    
    vector<ghostCell> ghostCells = data.ghostCells;
    
    for(int j=0; j<ghostCells.size(); j++)
    {
        o_latticeIndices.push_back(ghostCells[j].globalLatticeIdx);
        o_scalars.push_back(ghostCells[j].value);
    }
    
    vector<double> values = data.scalarValues;
    
    for(int j=0; j<values.size(); j++)
    {
        int glIdx = data.globalLatticeIndices[j];
        
        o_latticeIndices.push_back(glIdx);
        o_scalars.push_back(values[j]);
    }
    
    vector<double> i_scalars;
    vector<int> i_latticeIndices;
    
    for(int root=0; root<PAR_Size(); root++)
    {
        if(myRank == root)
        {
            int idxTot = o_latticeIndices.size();
            
            if(idxTot == 0)
            {
                o_latticeIndices.push_back(-1);
            }
            
            // Send the list of ghost lattice cell indices this rank has
            BroadcastIntVectorFromAny(o_latticeIndices, myRank, root);
            
            if(idxTot > 0)
            {
                // Send scalars
                BroadcastDoubleVectorFromAny(o_scalars, myRank, root);
            }
        }
        else    // Receive ghost cells
        {
            // Ghost Cell Global Lattice Indices
            BroadcastIntVectorFromAny(i_latticeIndices, myRank, root);
            
            if(i_latticeIndices[0] != -1)
            {
                // Receive scalars
                BroadcastDoubleVectorFromAny(i_scalars, myRank, root);
                
                // Update scalar values
                for(int j=0; j<i_latticeIndices.size(); j++)
                {
                    int localIdx = GetLocalIndex(data, i_latticeIndices[j]);
                    if(localIdx >= 0)
                    {
                        data.scalarValues[localIdx] += i_scalars[j];
                    }
                    else
                    {
                        int gcIdx = GetGhostCellIndex(data.ghostCells, i_latticeIndices[j]);
                        if(gcIdx >= 0)
                        {
                            data.ghostCells[gcIdx].value += i_scalars[j];
                        }
                    }
                }
            }
        }
    }
}

// ****************************************************************************
// Method:  avtSPHResampleFilter::ExtendBoundsIfNeeded
//
// Purpose:     Extend the dataset bounds, if necessary, to match the user supplied
//              min/max ranges.
//
// Arguments:
//      bounds  the dataset's current bounds
//
// Programmer:  Kevin Griffin
// Creation:    Wed Aug 17 16:17:16 PDT 2016
//
//  Modifications:
//
// ****************************************************************************
void
avtSPHResampleFilter::ExtendBoundsIfNeeded(double *const bounds)
{
    int myRank = PAR_Rank();
    bool skip = false;
    bool needMinX = true;
    bool needMaxX = true;
    bool needMinY = true;
    bool needMaxY = true;
    bool needMinZ = (nDim == 3) ? true : false;
    bool needMaxZ = (nDim == 3) ? true : false;
    
    if(bounds[0] == 0 && bounds[1] == 0 && bounds[2] == 0 && bounds[3] == 0 && bounds[4] == 0 && bounds[5] == 0)
    {
        skip = true;
    }
        
    double temp[6];
    temp[0] = bounds[0];
    temp[1] = bounds[1];
    temp[2] = bounds[2];
    temp[3] = bounds[3];
    temp[4] = bounds[4];
    temp[5] = bounds[5];
    
    for(int root=0; root<PAR_Size(); root++)
    {
        if(myRank == root)
        {
            // Send this node's dataset bounds
            BroadcastDoubleArrayFromAny(bounds, 6, myRank);
        }
        else
        {
            double recvBounds[6];
            
            // Receive dataset bounds from root
            BroadcastDoubleArrayFromAny(recvBounds, 6, root);
            
            if(!skip)
            {
                // Extend X dimensions if needed
                if(needMinX)
                {
                    if(recvBounds[0] < bounds[0])
                    {
                        if(bounds[2] >= recvBounds[2] && bounds[3] <= recvBounds[3])
                        {
                            if(nDim == 3)
                            {
                                if(bounds[4] >= recvBounds[4] && bounds[5] <= recvBounds[5])
                                {
                                    bounds[0] = temp[0];
                                    needMinX = false;
                                }
                            }
                            else
                            {
                                bounds[0] = temp[0];
                                needMinX = false;
                            }
                        }
                        else
                        {
                            bounds[0] = atts.GetMinX();
                        }
                    }
                }
                
                if(needMaxX)
                {
                    if(recvBounds[1] > bounds[1])
                    {
                        if(bounds[2] >= recvBounds[2] && bounds[3] <= recvBounds[3])
                        {
                            if(nDim == 3)
                            {
                                if(bounds[4] >= recvBounds[4] && bounds[5] <= recvBounds[5])
                                {
                                    bounds[1] = temp[1];
                                    needMaxX = false;
                                    
                                }
                            }
                            else
                            {
                                bounds[1] = temp[1];
                                needMaxX = false;
                            }
                        }
                        else
                        {
                            bounds[1] = atts.GetMaxX();
                        }
                    }
                }
                
                // Extend Y dimensions if needed
                if(needMinY)
                {
                    if(recvBounds[2] < bounds[2])
                    {
                        if(bounds[0] >= recvBounds[0] && bounds[1] <= recvBounds[1])
                        {
                            if(nDim == 3)
                            {
                                if(bounds[4] >= recvBounds[4] && bounds[5] <= recvBounds[5])
                                {
                                    bounds[2] = temp[2];
                                    needMinY = false;
                                }
                            }
                            else
                            {
                                bounds[2] = temp[2];
                                needMinY = false;
                            }
                        }
                        else
                        {
                            bounds[2] = atts.GetMinY();
                        }
                    }
                }
                
                if(needMaxY)
                {
                    if(recvBounds[3] > bounds[3])
                    {
                        if(bounds[0] >= recvBounds[0] && bounds[1] <= recvBounds[1])
                        {
                            if(nDim == 3)
                            {
                                if(bounds[4] >= recvBounds[4] && bounds[5] <= recvBounds[5])
                                {
                                    bounds[3] = temp[3];
                                    needMaxY = false;
                                }
                            }
                            else
                            {
                                bounds[3] = temp[3];
                                needMaxY = false;
                            }
                        }
                        else
                        {
                            bounds[3] = atts.GetMaxY();
                        }
                    }
                }
                
                // Extend Z dimensions if needed
                if(needMinZ)
                {
                    if(recvBounds[4] < bounds[4])
                    {
                        if(bounds[0] >= recvBounds[0] && bounds[1] <= recvBounds[1])
                        {
                            if(bounds[2] >= recvBounds[2] && bounds[3] <= recvBounds[3])
                            {
                                bounds[4] = temp[4];
                                needMinZ = false;
                            }
                        }
                        else
                        {
                            bounds[4] = atts.GetMinZ();
                        }
                    }
                }
                
                if(needMaxZ)
                {
                    if(recvBounds[5] > bounds[5])
                    {
                        if(bounds[0] >= recvBounds[0] && bounds[1] <= recvBounds[1])
                        {
                            if(bounds[2] >= recvBounds[2] && bounds[3] <= recvBounds[3])
                            {
                                bounds[5] = temp[5];
                                needMaxZ = false;
                            }
                        }
                        else
                        {
                            bounds[5] = atts.GetMaxZ();
                        }
                    }
                }
            }
        }
    }
}

// ****************************************************************************
// Method:  avtSPHResampleFilter::SynchMoments
//
// Purpose: Perform a sum reduction of the moments across all ranks.
//
// Arguments:
//      data  The SPH data containing this rank's moment values
//
// Programmer:  Kevin Griffin
// Creation:    Thu Jun 16 19:07:57 PDT 2016
//
//  Modifications:
//
// ****************************************************************************
void
avtSPHResampleFilter::SynchMoments(sphData &data)
{
    int myRank = PAR_Rank();
    
    vector<double> o_m0;
    vector<vector<double> > o_m1;
    vector<tensorStruct> o_m2;
    vector<int> o_latticeIndices;
    
    // Collect all the data that need to be synchronized
    vector<ghostCell> ghostCells = data.ghostCells;
    
    for(int j=0; j<ghostCells.size(); j++)
    {
        o_latticeIndices.push_back(ghostCells[j].globalLatticeIdx);
        o_m0.push_back(ghostCells[j].m0);
        o_m1.push_back(ghostCells[j].m1);
        o_m2.push_back(ghostCells[j].m2);
    }
    
    for(int j=0; j<data.totalCells; j++) {
        o_latticeIndices.push_back(data.globalLatticeIndices[j]);
        o_m0.push_back(data.m0[j]);
        o_m1.push_back(data.m1[j]);
        o_m2.push_back(data.m2[j]);
    }
    
    vector<double> i_m0;
    vector<vector<double> > i_m1;
    vector<tensorStruct> i_m2;
    vector<int> i_latticeIndices;
    
    for(int root=0; root<PAR_Size(); root++)
    {
        if(myRank == root)
        {
            int idxTot = o_latticeIndices.size();
           
            if(idxTot == 0) {
                o_latticeIndices.push_back(-1);
            }
            
            // Send the list of ghost lattice cell indices this rank has
            BroadcastIntVectorFromAny(o_latticeIndices, myRank, root);
            
            if(idxTot > 0)
            {
                // Send m0
                BroadcastDoubleVectorFromAny(o_m0, myRank, root);
                
                // Send m1
                double *m1Array = new double[nDim*idxTot];
                for(int i=0; i<idxTot; i++)
                {
                    for(int j=0; j<nDim; j++)
                    {
                        m1Array[i+nDim+j] = o_m1[i].at(j);
                    }
                }
                
                BroadcastDoubleArrayFromAny(m1Array, idxTot*nDim, root);
                delete [] m1Array;
                
                // Send m2
                BroadcastDoubleArrayFromAny(&o_m2[0].xx, idxTot*9, root);
            }
        }
        else    // Receive ghost cells
        {
            // Ghost Cell Global Lattice Indices
            BroadcastIntVectorFromAny(i_latticeIndices, myRank, root);
            int idxTot = i_latticeIndices.size();
            
            if(i_latticeIndices[0] != -1)
            {
                // Receive m0
                BroadcastDoubleVectorFromAny(i_m0, myRank, root);
                
                // Receive m1
                double *temp_m1 = new double[nDim*idxTot];
                i_m1.resize(idxTot);
                
                BroadcastDoubleArrayFromAny(temp_m1, idxTot*nDim, root);
                
                for(int i=0; i<idxTot; i++)
                {
                    vector<double> m1Val;
                    
                    for(int j=0; j<nDim; j++)
                    {
                        m1Val.push_back(temp_m1[i*nDim+j]);
                    }
                    
                    i_m1[i] = m1Val;
                }
                
                delete [] temp_m1;
                
                // Receive m2
                i_m2.resize(9*idxTot);
                BroadcastDoubleArrayFromAny(&i_m2[0].xx, 9*idxTot, root);
                
                // Update moment values
                for(int j=0; j<i_latticeIndices.size(); j++)
                {
                    int localIdx = GetLocalIndex(data, i_latticeIndices[j]);
                    if(localIdx >= 0)
                    {
                        data.m0[localIdx] += i_m0[j];
                        
                        if(nDim == 2)
                        {
                            data.m1[localIdx] = vectorSum<2>(data.m1[localIdx], i_m1[j]);
                            data.m2[localIdx] = tensorSum<2>(data.m2[localIdx], i_m2[j]);
                        }
                        else
                        {
                            data.m1[localIdx] = vectorSum<3>(data.m1[localIdx], i_m1[j]);
                            data.m2[localIdx] = tensorSum<3>(data.m2[localIdx], i_m2[j]);
                        }
                    }
                    else
                    {
                        int gcIdx = GetGhostCellIndex(data.ghostCells, i_latticeIndices[j]);
                        
                        if(gcIdx >= 0)
                        {
                            data.ghostCells[gcIdx].m0 += i_m0[j];
                            
                            if(nDim == 2)
                            {
                                data.ghostCells[gcIdx].m1 = vectorSum<2>(data.ghostCells[gcIdx].m1, i_m1[j]);
                                data.ghostCells[gcIdx].m2 = tensorSum<2>(data.ghostCells[gcIdx].m2, i_m2[j]);
                            }
                            else
                            {
                                data.ghostCells[gcIdx].m1 = vectorSum<3>(data.ghostCells[gcIdx].m1, i_m1[j]);
                                data.ghostCells[gcIdx].m2 = tensorSum<3>(data.ghostCells[gcIdx].m2, i_m2[j]);
                            }
                        }
                    }
                }
            }
        }
    } // End For
}

// ****************************************************************************
// Method:  avtSPHResampleFilter::GetDimensions
//
// Purpose: Generate the xyz dimensions for a dataset with bounds.
//
// Arguments:
//      dims    The output dimensions of the grid
//      min     The output minimum x,y, and z coordinates
//      max     The output maximum x,y, and z coordinates
//      steps   The input steps per x, y, and z
//      bounds  The input bounds of the dataset
//
//
// Programmer:  Kevin Griffin
// Creation:    Thu Jun 23 16:02:34 PDT 2016
//
//
// ****************************************************************************
void
avtSPHResampleFilter::GetDimensions(vector<int> &dims, vector<double> &min, vector<double> &max, const vector<double> &steps, const double *const bounds)
{
    // X Dims
    min[0] = (bounds[0] >= atts.GetMinX() && bounds[0] <= atts.GetMaxX()) ? bounds[0] : atts.GetMinX();
    max[0] = (bounds[1] <= atts.GetMaxX() && bounds[1] >= atts.GetMinX()) ? bounds[1] : atts.GetMaxX();
    dims[0] = int(ceil((max[0] - min[0]) / steps[0]));
    
    // Y Dims
    min[1] = (bounds[2] >= atts.GetMinY() && bounds[2] <= atts.GetMaxY()) ? bounds[2] : atts.GetMinY();
    max[1] = (bounds[3] <= atts.GetMaxY() && bounds[3] >= atts.GetMinY()) ? bounds[3] : atts.GetMaxY();
    dims[1] = int(ceil((max[1] - min[1]) / steps[1]));
    
    // Z Dims
    if(nDim == 3)
    {
        min[2] = (bounds[4] >= atts.GetMinZ() && bounds[4] <= atts.GetMaxZ()) ? bounds[4] : atts.GetMinZ();
        max[2] = (bounds[5] <= atts.GetMaxZ() && bounds[5] >= atts.GetMinZ()) ? bounds[5] : atts.GetMaxZ();
        dims[2] = int(ceil((max[2] - min[2]) / steps[2]));
        
        if (dims[0] <= 0 || dims[1] <= 0 || dims[2] < 0)
        {
            EXCEPTION1(VisItException, "The grid to SPHResample on is degenerate."
                       "Make sure that the number of samples in each direction "
                       "is positive.");
        }
        
        debug5 << "SPHResampling onto grid of dimensions: " << dims[0] << ", " << dims[1] << ", " << dims[2] << endl;
    }
    else
    {
        if (dims[0] <= 0 || dims[1] <= 0)
        {
            EXCEPTION1(VisItException, "The grid to SPHResample on is degenerate."
                       "Make sure that the number of samples in each direction "
                       "is positive.");
        }
        
        debug5 << "SPHResampling onto grid of dimensions: " << dims[0] << ", " << dims[1] << endl;
    }
}

// ****************************************************************************
// Method:  avtSPHResampleFilter::CreateOutputGrid
//
// Purpose: Map the sampled grid to the output grid.
//
// Arguments:
//      bounds  The grid bounds
//      steps The dimensions in each direction
//
//
// Programmer:  Cody Raskin
// Creation:    July 1, 2015
//
//  Modifications:
//  Kevin Griffin, Mon Jan 4 17:51:20 PST 2016
//  Changed return type to vtkRectilinear Grid and added a check for cell
//  centered data.
//
//  Kevin Griffin, Fri May 27 16:20:35 PDT 2016
//  Generalized the method by adding bounds and steps parameters.
//
// ****************************************************************************

vtkRectilinearGrid *
avtSPHResampleFilter::CreateOutputGrid(const double* const bounds, const vector<double> &steps)
{
    vtkRectilinearGrid *out_dset = vtkRectilinearGrid::New();
    vector<double> minXYZ(nDim, 0);
    vector<double> maxXYZ(nDim, 0);
    vector<int> dimVec(nDim, 0);

    GetDimensions(dimVec, minXYZ, maxXYZ, steps, bounds);
    
    int dims[3] = {0, 0, 1};
    dims[0] = dimVec[0] + 1;
    dims[1] = dimVec[1] + 1;
    
    // X Dims
    double xmin = minXYZ[0];
    
    // Y Dims
    double ymin = minXYZ[1];
  
    // Z Dims
    double zmin = 0;
    if(nDim == 3)
    {
        zmin = minXYZ[2];
        dims[2] = dimVec[2] + 1;
    }
    
    debug5 << "Output Grid - dim1: " << dims[0] << " dim2: " << dims[1] << " dims3: " << dims[2] << endl;
    debug5 << "Output Grid - step1: " << steps[0] << " step2: " << steps[1] << " step3: " << steps[2] << endl;
    
    out_dset->SetDimensions(dims);
    
    // create the xcoords
    vtkDoubleArray *out_xcoords = vtkDoubleArray::New();
    out_xcoords->SetNumberOfTuples(dims[0]);
    out_xcoords->SetNumberOfComponents(1);
    
    double dx = steps[0];
    
    for(int i=0;i<dims[0];i++)
    {
        double val = xmin + i * dx;
        val = val <= atts.GetMaxX() ? val : atts.GetMaxX();
        out_xcoords->SetTuple1(i,val);
    }
    
    out_dset->SetXCoordinates(out_xcoords);
    out_xcoords->Delete();
    
    // create the ycoords
    vtkDoubleArray *out_ycoords = vtkDoubleArray::New();
    out_ycoords->SetNumberOfTuples(dims[1]);
    out_ycoords->SetNumberOfComponents(1);

    double dy = steps[1];
    
    for(int i=0;i<dims[1];i++)
    {
        double val = ymin + i * dy;
        val = val <= atts.GetMaxY() ? val : atts.GetMaxY();
        out_ycoords->SetTuple1(i,val);
    }
    
    out_dset->SetYCoordinates(out_ycoords);
    out_ycoords->Delete();
    
    // create the zcoords
    vtkDoubleArray *out_zcoords = vtkDoubleArray::New();
    out_zcoords->SetNumberOfComponents(1);
    
    if(nDim == 3)
    {
        out_zcoords->SetNumberOfTuples(dims[2]);
        double dz   = steps[2];
        
        for(int i=0;i<dims[2];i++)
        {
            double val = zmin + i * dz;
            val = val <= atts.GetMaxZ() ? val : atts.GetMaxZ();
            out_zcoords->SetTuple1(i,val);
        }
    }
    else
    {
        out_zcoords->SetNumberOfTuples(1);
        out_zcoords->SetTuple1(0,zmin);
    }
    
    out_dset->SetZCoordinates(out_zcoords);
    out_zcoords->Delete();
    
    return out_dset;
}

// ****************************************************************************
//  Method: avtSPHResampleFilter::GetLocalIndex
//
//  Purpose:
//      Convert the global lattice index to the local lattice index. If the
//      global lattice index is out of range for the local lattice -1 is
//      is returned.
//
//  Arguments:
//      data        The SPH data (moments, correction coefficients, etc)
//      globalIdx   The global lattice index
//
//  Returns:        The local lattice index corresponding to the global lattice
//                  index, otherwise -1
//
//  Programmer: Kevin Griffin
//  Creation:   Thu May 26 16:20:35 PDT 2016
//
//  Modifications:
//
// ****************************************************************************
int
avtSPHResampleFilter::GetLocalIndex(const sphData &data, int globalIdx)
{
    if(nDim == 2)
    {
        int y = globalIdx / atts.GetXnum();
        int x = globalIdx - (y * atts.GetXnum());
        
        if((x >= data.indices[0] && x < data.indices[1]) && (y >= data.indices[2] && y < data.indices[3]))
        {
            int localx = x - data.indices[0];
            int localy = y - data.indices[2];
            
            int xNum = data.indices[1] - data.indices[0];
            
            int localIdx = localy * xNum + localx;
            return localIdx;
        }
    }
    else    // nDim == 3
    {
        int z = globalIdx / (atts.GetXnum() * atts.GetYnum());
        int tempIdx = globalIdx - (z * atts.GetXnum() * atts.GetYnum());
        int y =  tempIdx / atts.GetXnum();
        int x = tempIdx % atts.GetXnum();
        
        if((x >= data.indices[0] && x < data.indices[1]) && (y >= data.indices[2] && y < data.indices[3]) && (z >= data.indices[4] && z < data.indices[5]))
        {
            int localx = x - data.indices[0];
            int localy = y - data.indices[2];
            int localz = z - data.indices[4];
            
            int yNum = data.indices[3] - data.indices[2];
            int xNum = data.indices[1] - data.indices[0];
            
            int localIdx = localz * yNum * xNum + localy * xNum + localx;
            return localIdx;
        }
    }
    
    return -1;
}

// ****************************************************************************
//  Method: avtSPHResampleFilter::GetGhostCellIndex
//
//  Purpose:
//      Find the index of the ghost cell with global lattice index equal to
//      globalIdx.
//
//  Arguments:
//      ghostCells  The list of ghost cells
//      globalIdx   The global lattice index
//
//  Returns:       The index of the ghost cell in the list, otherwise -1
//
//  Programmer: Kevin Griffin
//  Creation:   Thu May 26 16:20:35 PDT 2016
//
//  Modifications:
//
// ****************************************************************************
int
avtSPHResampleFilter::GetGhostCellIndex(const vector<ghostCell> &ghostCells,  int globalIdx)
{
    for(int i=0; i<ghostCells.size(); i++)
    {
        if(ghostCells[i].globalLatticeIdx == globalIdx)
        {
            return i;
        }
    }
    
    return -1;
}

// ****************************************************************************
// Method:  avtSPHResampleFilter::CheckInputVariables
//
// Purpose: Check that the required variables exist and have the same centering.
//
// Arguments:
//      var1    resample variable
//      var2    support variable
//      var3    weight variable
//
//
// Programmer:  Kevin Griffin
// Creation:    TWed May 25 14:38:01 PDT 2016
//
//  Modifications:
//
// ****************************************************************************
void
avtSPHResampleFilter::CheckInputVariables(const vtkDataArray *const var1,
                                          const vtkDataArray *const var2,
                                          const vtkDataArray *const var3)
{
    std::string errMsg;
    
    if(var1 == NULL) {
        errMsg.append(resampleVarName);
        errMsg.append(" must exist and have the same centering");
        EXCEPTION1(InvalidVariableException, errMsg);
    }
    
    if(var2 == NULL) {
        errMsg.clear();
        errMsg.append(supportVarName);
        errMsg.append(" must exist and have the same centering");
        EXCEPTION1(InvalidVariableException, errMsg);
    }
    
    if(var3 == NULL) {
        errMsg.clear();
        errMsg.append(weightVarName);
        errMsg.append(" must exist and have the same centering");
        EXCEPTION1(InvalidVariableException, errMsg);
    }
}

// ****************************************************************************
//  Method: avtSPHResampleFilter::ModifyContract
//
//  Purpose:
//      Creates a contract the removes the operator-created-expression.
//
//  Programmer: harrison37 -- generated by xml2avt
//  Creation:   Fri Dec 5 13:50:31 PST 2014
//
//  Modifications:
//  Kevin Griffin, Mon Jan 4 17:51:20 PST 2016
//  Removed the avtResampleSelection code.
//
// ****************************************************************************

avtContract_p
avtSPHResampleFilter::ModifyContract(avtContract_p in_contract)
{
    avtContract_p rv = new avtContract(in_contract);
    rv->GetDataRequest()->SetDesiredGhostDataType(NO_GHOST_DATA);   // GHOST_ZONE_DATA
    rv->NoStreaming();
    
    resampleVarName = std::string(rv->GetDataRequest()->GetVariable());
    supportVarName  = atts.GetTensorSupportVariable();
    weightVarName   = atts.GetWeightVariable();
    
    if( supportVarName == "default")
    {
        EXCEPTION1(InvalidVariableException,"default");
    }
    
    if( weightVarName == "default")
    {
        weightVarName = "mass";
    }
   
    debug5 << "Resample Variable Name: " << resampleVarName << endl;
    debug5 << "Tensor Support Variable Name: " << supportVarName << endl;
    debug5 << "Weight Variable Name: " << weightVarName << endl;
    
    // make sure to request the tensor support var
    rv->GetDataRequest()->AddSecondaryVariable(supportVarName.c_str());
    rv->GetDataRequest()->AddSecondaryVariable(weightVarName.c_str());
    
    if (in_contract->GetDataRequest()->MayRequireZones() ||
        in_contract->GetDataRequest()->MayRequireNodes())
    {
        avtDataAttributes &data = GetInput()->GetInfo().GetAttributes();
        
        if (data.ValidActiveVariable())
        {
            if (data.GetCentering() == AVT_NODECENT)
            {
                rv->GetDataRequest()->TurnNodeNumbersOn();
            }
            else if (data.GetCentering() == AVT_ZONECENT)
            {
                rv->GetDataRequest()->TurnZoneNumbersOn();
            }
        }
        else
        {
            // canot determine variable centering, so turn on both
            // node numbers and zone numbers.
            rv->GetDataRequest()->TurnNodeNumbersOn();
            rv->GetDataRequest()->TurnZoneNumbersOn();
        }
    }
    
    return rv;
}

// ****************************************************************************
//  Method: avtSPHResampleFilter::UpdateDataObjectInfo
//
//  Purpose:
//      Tells output that we have a new variable.
//
//  Programmer: harrison37 -- generated by xml2avt
//  Creation:   Fri Dec 5 13:50:31 PST 2014
//
//  Modifications:
//  Kevin Griffin, Mon Jan 4 17:51:20 PST 2016
//  Added metadata, desired spatial/data extents, and commented out
//  InvalidateNodes.
//
// ****************************************************************************

void
avtSPHResampleFilter::UpdateDataObjectInfo(void)
{
    avtDataAttributes &inAtts   = GetInput()->GetInfo().GetAttributes();
    avtDataAttributes &outAtts = GetOutput()->GetInfo().GetAttributes();
    int spatialDim = inAtts.GetSpatialDimension();
    
    GetOutput()->GetInfo().GetValidity().InvalidateZones();
    outAtts.SetTopologicalDimension(spatialDim);
    
    double bounds[6] = {0, 0, 0, 0, 0, 0};
    bounds[0] = atts.GetMinX();
    bounds[1] = atts.GetMaxX();
    bounds[2] = atts.GetMinY();
    bounds[3] = atts.GetMaxY();
    bounds[4] = atts.GetMinZ();
    bounds[5] = spatialDim == 3 ? atts.GetMaxZ() : atts.GetMinZ();
    outAtts.GetDesiredSpatialExtents()->Set(bounds);
    
    if (!resampleVarName.empty() && outAtts.ValidActiveVariable())
    {
        if (outAtts.GetVariableName() != resampleVarName)
            outAtts.SetActiveVariable(resampleVarName.c_str());
        
        double range[2];
        GetDataExtents(range, resampleVarName.c_str());
        outAtts.GetDesiredDataExtents()->Set(range);
    }
    
    outAtts.RemoveVariable(supportVarName);
    outAtts.RemoveVariable(weightVarName);
    
    char params[200];
    snprintf(params, 200, "nx=%d ny=%d nz=%d", atts.GetXnum(), atts.GetYnum(), (spatialDim == 3 ? atts.GetZnum() : 0));
    outAtts.AddFilterMetaData("SPH Resample", params);
}
